function test_matrices()
% TEST_MATRICES LU分解测试矩阵生成器
% 
% 生成各种类型的测试矩阵用于LU分解算法测试

fprintf('=== LU分解测试矩阵生成器 ===\n\n');

%% 1. 随机矩阵
fprintf('1. 随机矩阵\n');

sizes = [10, 50, 100];
fprintf('生成不同尺寸的随机矩阵:\n');
fprintf('尺寸    条件数        行列式        特征\n');
fprintf('----    ------        ------        ----\n');

random_matrices = {};
for i = 1:length(sizes)
    n = sizes(i);
    rng(42 + i);  % 固定随机种子
    A = randn(n, n);
    
    cond_num = cond(A);
    det_val = det(A);
    
    fprintf('%4d    %8.2e    %10.2e    一般随机矩阵\n', n, cond_num, det_val);
    
    random_matrices{i} = A;
end

%% 2. 病态矩阵
fprintf('\n2. 病态矩阵\n');

fprintf('生成高条件数的病态矩阵:\n');
fprintf('类型              尺寸    条件数        特征\n');
fprintf('----------------  ----    ------        ----\n');

ill_conditioned_matrices = {};

% Hilbert矩阵
for n = [5, 8, 10]
    H = hilb(n);
    cond_H = cond(H);
    fprintf('Hilbert矩阵       %4d    %8.2e    经典病态矩阵\n', n, cond_H);
    ill_conditioned_matrices{end+1} = H;
end

% Vandermonde矩阵
for n = [8, 10, 12]
    x = linspace(0, 1, n);
    V = vander(x);
    cond_V = cond(V);
    fprintf('Vandermonde矩阵   %4d    %8.2e    多项式基矩阵\n', n, cond_V);
    ill_conditioned_matrices{end+1} = V;
end

% 人工构造的病态矩阵
for n = [10, 20]
    [U, ~, V] = svd(randn(n));
    s = logspace(0, -12, n);  % 指数衰减的奇异值
    A_ill = U * diag(s) * V';
    cond_A = cond(A_ill);
    fprintf('人工病态矩阵      %4d    %8.2e    指数衰减奇异值\n', n, cond_A);
    ill_conditioned_matrices{end+1} = A_ill;
end

%% 3. 特殊结构矩阵
fprintf('\n3. 特殊结构矩阵\n');

special_matrices = {};

% 三对角矩阵
fprintf('三对角矩阵:\n');
for n = [10, 20, 50]
    A_tri = diag(2*ones(n,1)) + diag(-ones(n-1,1),1) + diag(-ones(n-1,1),-1);
    cond_tri = cond(A_tri);
    fprintf('  尺寸%d: 条件数 = %.2e\n', n, cond_tri);
    special_matrices{end+1} = A_tri;
end

% 带状矩阵
fprintf('带状矩阵 (带宽=3):\n');
for n = [15, 30]
    A_band = diag(4*ones(n,1)) + diag(-ones(n-1,1),1) + diag(-ones(n-1,1),-1) + ...
             diag(0.5*ones(n-2,1),2) + diag(0.5*ones(n-2,1),-2);
    cond_band = cond(A_band);
    fprintf('  尺寸%d: 条件数 = %.2e\n', n, cond_band);
    special_matrices{end+1} = A_band;
end

% 循环矩阵
fprintf('循环矩阵:\n');
for n = [8, 16]
    c = [4, -1, 0, 0, 0, 0, 0, -1];  % 第一行
    if n > 8
        c = [c, zeros(1, n-8)];
    else
        c = c(1:n);
    end
    A_circ = gallery('circul', c);
    cond_circ = cond(A_circ);
    fprintf('  尺寸%d: 条件数 = %.2e\n', n, cond_circ);
    special_matrices{end+1} = A_circ;
end

%% 4. 稀疏矩阵
fprintf('\n4. 稀疏矩阵\n');

sparse_matrices = {};

fprintf('稀疏矩阵生成:\n');
fprintf('类型        尺寸    稀疏度    非零元素    条件数\n');
fprintf('--------    ----    ------    --------    ------\n');

% 随机稀疏矩阵
for n = [50, 100]
    density = 0.1;  % 10%非零元素
    A_sparse = sprandn(n, n, density) + speye(n);  % 加对角元素确保非奇异
    A_sparse_full = full(A_sparse);
    
    nnz_count = nnz(A_sparse);
    sparsity = 1 - nnz_count / n^2;
    cond_sparse = cond(A_sparse_full);
    
    fprintf('随机稀疏    %4d    %6.3f    %8d    %8.2e\n', ...
            n, sparsity, nnz_count, cond_sparse);
    sparse_matrices{end+1} = A_sparse_full;
end

% 五点差分矩阵（来自偏微分方程）
for grid_size = [10, 15]
    n = grid_size^2;
    A_pde = gallery('poisson', grid_size);
    A_pde_full = full(A_pde);
    
    nnz_count = nnz(A_pde);
    sparsity = 1 - nnz_count / n^2;
    cond_pde = cond(A_pde_full);
    
    fprintf('五点差分    %4d    %6.3f    %8d    %8.2e\n', ...
            n, sparsity, nnz_count, cond_pde);
    sparse_matrices{end+1} = A_pde_full;
end

%% 5. 正定矩阵
fprintf('\n5. 正定矩阵\n');

positive_definite_matrices = {};

fprintf('正定矩阵生成:\n');
fprintf('类型            尺寸    条件数        最小特征值\n');
fprintf('------------    ----    ------        ----------\n');

% 对称正定随机矩阵
for n = [10, 20]
    rng(100 + n);
    A_temp = randn(n, n);
    A_pd = A_temp' * A_temp + eye(n);  % 确保正定
    
    cond_pd = cond(A_pd);
    min_eig = min(eig(A_pd));
    
    fprintf('随机正定        %4d    %8.2e    %10.4f\n', n, cond_pd, min_eig);
    positive_definite_matrices{end+1} = A_pd;
end

% 相关矩阵
for n = [8, 15]
    % 生成相关矩阵
    rng(200 + n);
    X = randn(n, 2*n);  % 样本矩阵
    R = corrcoef(X');   % 相关矩阵
    
    cond_corr = cond(R);
    min_eig_corr = min(eig(R));
    
    fprintf('相关矩阵        %4d    %8.2e    %10.4f\n', n, cond_corr, min_eig_corr);
    positive_definite_matrices{end+1} = R;
end

%% 6. 近奇异矩阵
fprintf('\n6. 近奇异矩阵\n');

near_singular_matrices = {};

fprintf('近奇异矩阵生成:\n');
fprintf('类型              尺寸    条件数        最小奇异值\n');
fprintf('----------------  ----    ------        ----------\n');

% 秩亏矩阵的扰动
for n = [10, 15]
    rng(300 + n);
    % 构造秩亏矩阵
    A_rank_def = randn(n, n-2) * randn(n-2, n);
    % 添加小扰动
    epsilon = 1e-10;
    A_near_sing = A_rank_def + epsilon * randn(n, n);
    
    cond_ns = cond(A_near_sing);
    min_sv = min(svd(A_near_sing));
    
    fprintf('扰动秩亏矩阵    %4d    %8.2e    %10.2e\n', n, cond_ns, min_sv);
    near_singular_matrices{end+1} = A_near_sing;
end

% 几乎相关的行
for n = [8, 12]
    rng(400 + n);
    A_corr_rows = randn(n, n);
    % 使最后一行几乎等于前面行的线性组合
    coeffs = randn(n-1, 1);
    coeffs = coeffs / sum(abs(coeffs));  % 归一化
    A_corr_rows(end, :) = coeffs' * A_corr_rows(1:end-1, :) + 1e-12 * randn(1, n);
    
    cond_cr = cond(A_corr_rows);
    min_sv_cr = min(svd(A_corr_rows));
    
    fprintf('相关行矩阵      %4d    %8.2e    %10.2e\n', n, cond_cr, min_sv_cr);
    near_singular_matrices{end+1} = A_corr_rows;
end

%% 7. 实际应用矩阵
fprintf('\n7. 实际应用矩阵\n');

application_matrices = {};

% 离散化拉普拉斯算子
fprintf('应用矩阵生成:\n');
fprintf('应用领域        尺寸    条件数        特征\n');
fprintf('------------    ----    ------        ----\n');

for grid_size = [8, 12]
    n = grid_size^2;
    % 二维拉普拉斯算子的五点差分格式
    A_laplace = gallery('poisson', grid_size);
    A_laplace_full = full(A_laplace);
    
    cond_lap = cond(A_laplace_full);
    fprintf('偏微分方程      %4d    %8.2e    拉普拉斯算子\n', n, cond_lap);
    application_matrices{end+1} = A_laplace_full;
end

% 马尔可夫链转移矩阵
for n = [10, 15]
    rng(500 + n);
    % 生成随机转移矩阵
    P = rand(n, n);
    P = P ./ sum(P, 2);  % 行归一化
    % 转换为线性系统矩阵
    A_markov = eye(n) - P';
    A_markov(end, :) = ones(1, n);  % 添加归一化约束
    
    cond_markov = cond(A_markov);
    fprintf('马尔可夫链      %4d    %8.2e    转移矩阵\n', n, cond_markov);
    application_matrices{end+1} = A_markov;
end

%% 8. 测试矩阵性能评估
fprintf('\n8. 测试矩阵性能评估\n');

% 选择代表性矩阵进行LU分解测试
test_cases = {
    {'随机矩阵', random_matrices{2}},
    {'Hilbert矩阵', hilb(8)},
    {'三对角矩阵', special_matrices{2}},
    {'稀疏矩阵', sparse_matrices{1}},
    {'正定矩阵', positive_definite_matrices{1}},
    {'近奇异矩阵', near_singular_matrices{1}}
};

fprintf('LU分解性能测试:\n');
fprintf('矩阵类型        尺寸    LU时间(s)    误差        选主元必要性\n');
fprintf('------------    ----    ---------    ----        ------------\n');

for i = 1:length(test_cases)
    matrix_name = test_cases{i}{1};
    A = test_cases{i}{2};
    n = size(A, 1);
    
    % 不选主元的LU分解
    try
        tic;
        [L_no_pivot, U_no_pivot] = lu_no_pivot(A);
        time_no_pivot = toc;
        error_no_pivot = norm(A - L_no_pivot * U_no_pivot, 'fro') / norm(A, 'fro');
        pivot_needed = '否';
    catch
        time_no_pivot = inf;
        error_no_pivot = inf;
        pivot_needed = '是';
    end
    
    % 选主元的LU分解
    tic;
    [L_pivot, U_pivot, P_pivot] = lu(A);
    time_pivot = toc;
    error_pivot = norm(A - P_pivot' * L_pivot * U_pivot, 'fro') / norm(A, 'fro');
    
    % 选择更好的结果显示
    if error_no_pivot < 1e-10
        display_time = time_no_pivot;
        display_error = error_no_pivot;
    else
        display_time = time_pivot;
        display_error = error_pivot;
        pivot_needed = '是';
    end
    
    fprintf('%-12s    %4d    %9.6f    %8.2e    %s\n', ...
            matrix_name, n, display_time, display_error, pivot_needed);
end

%% 9. 矩阵生成建议
fprintf('\n9. 测试矩阵选择建议\n');

fprintf('不同测试目的的矩阵选择:\n\n');

fprintf('算法正确性测试:\n');
fprintf('• 小尺寸随机矩阵 (5×5 到 10×10)\n');
fprintf('• 已知解析解的矩阵\n');
fprintf('• 单位矩阵和对角矩阵\n\n');

fprintf('数值稳定性测试:\n');
fprintf('• Hilbert矩阵 (高条件数)\n');
fprintf('• Vandermonde矩阵\n');
fprintf('• 近奇异矩阵\n\n');

fprintf('性能测试:\n');
fprintf('• 大尺寸随机矩阵 (100×100 以上)\n');
fprintf('• 稀疏矩阵\n');
fprintf('• 带状矩阵\n\n');

fprintf('特殊结构测试:\n');
fprintf('• 正定矩阵 (Cholesky分解)\n');
fprintf('• 三对角矩阵 (专用算法)\n');
fprintf('• 循环矩阵\n\n');

fprintf('实际应用测试:\n');
fprintf('• 偏微分方程离散化矩阵\n');
fprintf('• 最小二乘问题的法方程\n');
fprintf('• 马尔可夫链转移矩阵\n\n');

%% 10. 保存测试矩阵
fprintf('\n10. 测试矩阵保存\n');

% 创建测试矩阵集合
test_matrix_collection = struct();
test_matrix_collection.random = random_matrices;
test_matrix_collection.ill_conditioned = ill_conditioned_matrices;
test_matrix_collection.special = special_matrices;
test_matrix_collection.sparse = sparse_matrices;
test_matrix_collection.positive_definite = positive_definite_matrices;
test_matrix_collection.near_singular = near_singular_matrices;
test_matrix_collection.application = application_matrices;

fprintf('测试矩阵集合已创建，包含:\n');
fprintf('• %d个随机矩阵\n', length(random_matrices));
fprintf('• %d个病态矩阵\n', length(ill_conditioned_matrices));
fprintf('• %d个特殊结构矩阵\n', length(special_matrices));
fprintf('• %d个稀疏矩阵\n', length(sparse_matrices));
fprintf('• %d个正定矩阵\n', length(positive_definite_matrices));
fprintf('• %d个近奇异矩阵\n', length(near_singular_matrices));
fprintf('• %d个应用矩阵\n', length(application_matrices));

fprintf('\n使用方法:\n');
fprintf('matrices = test_matrices();\n');
fprintf('A = matrices.random{1};  %% 获取第一个随机矩阵\n');
fprintf('[L, U, P] = lu(A);       %% 进行LU分解\n\n');

% 返回矩阵集合
if nargout > 0
    varargout{1} = test_matrix_collection;
end

end

% 辅助函数：不选主元的LU分解
function [L, U] = lu_no_pivot(A)
n = size(A, 1);
L = eye(n);
U = A;

for k = 1:n-1
    if abs(U(k, k)) < eps
        error('需要选主元');
    end
    
    for i = k+1:n
        L(i, k) = U(i, k) / U(k, k);
        U(i, k:n) = U(i, k:n) - L(i, k) * U(k, k:n);
    end
end
end